Ссылка на практикум
from IPython.display import Image
import os, sys
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd
from IPython.display import Image
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
Информация о моделировании:
* Силовое поле используемое при построении: amber99sb * Заряд системы: 0 (-10+10) * Размер и форму ячейки: прямоугольник 5.092 4.906 5.462 (nm) * Минимизация энергии: - Алгоритм минимизации энергии: l-bfgs (steepest descent algoritm) - Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: * Модель, которой описывался растворитель: * Утряска растворителя: - Параметр который обуславливает неподвижность биополимера: - Число шагов: 1000 - Длина шага: 0.0002 ps - Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: - Алгоритм термостата: V-rescale - Алгоритм баростата: Berendsen * Основной расчёт МД - Время моделирования: 10 ns - Количество процессоров: 1 (?) - Эффективность масштабирования: - Длина траектории: - Число шагов: 5000000 - Длина шага: 0.002 ps - Алгоритм интегратора: md - Алгоритмы расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: pme и Cut-off - Алгоритм термостата: V-rescale - Алгоритм баростата: Berendsen
Визуальный анализ движений молекул (из консоли, выбирая группу №1):
%%bash
trjconv -f dna_md.xtc -s dna_md.tpr -o dna_pbc_1.pdb -skip 20 -pbc mol
%%bash
# с центрированием
trjconv -f dna_md.xtc -s dna_md.tpr -o dna_fit_2.pdb -skip 20 -fit rot+trans
%%bash
pymol dna_pbc_1.pdb
Кадр №1:
Image(filename="A-form.png")
Кадр №51 (последний):
Image(filename="B-form.png")
Сложно сказать, когда происходит переход форм ДНК, но скорее всего это просходит между кадрами 36-37.
Определение средне-квадратичное отколнение в ходе моделирования¶
Сначала расчитаем отклонение в ходе всей симуляции относительно стартовой структуры:
%%bash
g_rms -f dna_md.xtc -s dna_md.tpr -o rms_1
%%bash
sed -i 's/@/#/' rms_1.xvg
a = np.loadtxt("rms_1.xvg",comments='#')
x = a[:,0]
y = a[:,1]
#plt.rcParams["figure.figsize"] = [15,8]
plt.scatter(x, y)
plt.title('RMSD')
plt.suptitle('DNA after lsq fit to DNA')
plt.xlabel('Time (ps)')
plt.ylabel('RMSD (nm)')
plt.show()
Резкие скачки RMSD скорее всего соответствуют случаям, когда цепи ДНК разлетаются друг от друга (в разные ячейки).
В диапазоне 4000-6000 ps скорее всего произошёл переход ДНК из А в В форму, т.к. после разрыва (в базовой линии) базальный уровень колебаний (RMSD) стал немного больше, чем был, что может говорить о том, что цепи ДНК немного разошлись и получили большую свободу для колебаний.
Теперь - относительно каждой предыдущей структуры на растоянии 400 кадров:
%%bash
g_rms -f dna_md.xtc -s dna_md.tpr -o rms_2 -prev 400
%%bash
sed -i 's/@/#/' rms_2.xvg
a = np.loadtxt("rms_2.xvg",comments='#')
x = a[:,0]
y = a[:,1]
plt.scatter(x, y)
plt.title('RMSD')
plt.suptitle('RMSD with frame 4000 Time (ps) ago')
plt.xlabel('Time (ps)')
plt.ylabel('RMSD (nm)')
plt.show()
Анализ поверхности¶
Определим изменение гидрофобной и гидрофильной поверхности в ходе конформационного перехода:
%%bash
g_sas -f dna_md.xtc -s dna_md.tpr -o sas_dna.xvg
%%bash
sed -i 's/@/#/' sas_dna.xvg
a = np.loadtxt("sas_dna.xvg",comments='#')
x = a[:,0]
y1 = a[:,1]
y2 = a[:,2]
y3 = a[:,3]
y4 = a[:,4]
plt.scatter(x, y1, c='yellow')
plt.scatter(x, y2, c='blue')
#plt.scatter(x, y3, c='red') #Total
#plt.scatter(x, y4, c='grey') #D Gsolv
plt.legend(['Hydrophobic', "Hydrophilic"], loc=7)
plt.title("Solvent Accessible Surface")
plt.xlabel('Time (ps)')
plt.ylabel('Area (nm\S2\N)')
plt.show()
Видно, что к концу моделирования площадь гидрофобной поверхности немного уменьшилась. Резких скачков не происходило. Возможно, это способствовало переходу ДНК из А-формы в В-форму.
Анализ водородных связей¶
Исследуем связи между ДНК и ДНК, то есть водородные связи между цепями ДНК:
%%bash
g_hbond -f dna_md.xtc -s dna_md.tpr -num hbond_dna
%%bash
sed -i 's/@/#/' hbond_dna.xvg
a = np.loadtxt("hbond_dna.xvg",comments='#')
x = a[:,0]
y1 = a[:,1]
y2 = a[:,2]
plt.scatter(x, y1, c='blue')
plt.scatter(x, y2, c='magenta')
z = np.polyfit(x, y2, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'r-', linewidth=2)
plt.legend(['Trendline', 'Hydrogen bonds', "Pairs within 0.35 nm"], bbox_to_anchor = (1.6, 1))
plt.title("Hydrogen Bonds")
plt.xlabel('Time (ps)')
plt.ylabel('Number')
plt.show()
Интересно, что количество водородных связей между цепями ДНК уменьшилось. К тому же, это произошло из-за того, что увеличилось расстояние между атомами (trendline). Следовательно, цепи ДНК как бы раздвинулись, чтобы увеличить количество взаимодействий с окружающей водой (см. ниже). На самом деле, этого следовало ожидать, т.к. А-форма ДНК преобладает в сухой ДНК (и в РНК), а в окружении воды ДНК переходит в В-форму.
Изучим количество водородных связей ДНК-Вода:
%%bash
g_hbond -f dna_md.xtc -s dna_md.tpr -num hbond_dna_sol
%%bash
sed -i 's/@/#/' hbond_dna_sol.xvg
a = np.loadtxt("hbond_dna_sol.xvg",comments='#')
x = a[:,0]
y1 = a[:,1]
y2 = a[:,2]
plt.scatter(x, y1, c='blue')
plt.scatter(x, y2, c='magenta')
z = np.polyfit(x, y1, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'r-', linewidth=2)
plt.legend(['Trendline', 'Hydrogen bonds', "Pairs within 0.35 nm"], bbox_to_anchor = (1.6, 1))
plt.title("Hydrogen Bonds")
plt.xlabel('Time (ps)')
plt.ylabel('Number')
plt.show()
Предположение, сделанное выше, действительно подтвердилось.